MODULE sefmm_cheb_Fij

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: cheb_init_Fij,   &
             cheb_free_Fij,   &
             cheb_get_Fij

   TYPE(sefmm_Fij_array), SAVE :: Fvec(0:1,128,128)
   
   INTEGER(INTK), SAVE :: quad_pts = 300
   REAL(REALK),   SAVE :: T_VALUE_THR = 50.0d0
!      INTEGER(INTK), PARAMETER :: quad_pts = 1000
!      REAL(REALK),   PARAMETER :: T_VALUE_THR = 1000.0d0

CONTAINS

!------------------------------------------------------------------------------

   SUBROUTINE cheb_init_Fij(scheme)

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN) :: scheme

      INTEGER(INTK) :: i,j,k

      quad_pts =  scheme%Fij_nquad
      T_VALUE_THR = scheme%Fij_t_thresh

      DO i = 1, SIZE(Fvec,3)
         DO j = 1, SIZE(Fvec,2)
            DO k = 0, 1
               IF (ASSOCIATED(Fvec(k,j,i)%ij)) CALL fmm_quit('cheb_init_Fij')
               NULLIFY(Fvec(k,j,i)%ij)
            END DO
         END DO
      END DO

   END SUBROUTINE cheb_init_Fij

!------------------------------------------------------------------------------

   SUBROUTINE cheb_free_Fij

      IMPLICIT NONE
      INTEGER(INTK) :: i,j,k

      DO i = 1, SIZE(Fvec,3)
         DO j = 1, SIZE(Fvec,2)
            DO k = 0, 1
               IF (ASSOCIATED(Fvec(k,j,i)%ij)) DEALLOCATE (Fvec(k,j,i)%ij)
               NULLIFY(Fvec(k,j,i)%ij)
            END DO
         END DO
      END DO

   END SUBROUTINE cheb_free_Fij

!------------------------------------------------------------------------------

   SUBROUTINE cheb_get_Fij(nnbox,nrho,npot,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK),        INTENT(IN)  :: nnbox, nrho, npot
      REAL(REALK),          INTENT(IN)  :: t(:)
      REAL(REALK), POINTER, INTENT(OUT) :: Fij(:,:,:)

      LOGICAL :: fail
      INTEGER(INTK) :: order

      order = npot/(1+2*nnbox)

      ! If Fij is not already built, allocate memory and build it...
      IF (.NOT. ASSOCIATED(Fvec(nnbox,nrho,order)%ij)) THEN
!         WRITE(LUPRI,*) 'New F-vector, order(rho,V) = ', nrho, order
!         WRITE(LUPRI,*) 'Allocating F-vector...',          & 
!            (nrho*npot*SIZE(t))*8d0/(1024*1024), 'MB'
         ALLOCATE( Fvec(nnbox,nrho,order)%ij(nrho,npot,SIZE(t)) )
         CALL cheb_build_Fij(nnbox,nrho,npot,t,Fvec(nnbox,nrho,order)%ij)
      END IF

      Fij => Fvec(nnbox,nrho,order)%ij(:,:,:)

      ! Error checks
      fail = .FALSE.
      !print *, nnbox, nrho, order, npot
      IF (SIZE(Fij,1) /= nrho) fail = .TRUE.
      IF (SIZE(Fij,2) /= npot) fail = .TRUE.
      IF (SIZE(Fij,3) /= SIZE(t)) fail = .TRUE.
      IF (fail) CALL fmm_quit('error in Fij build')

   END SUBROUTINE cheb_get_Fij

!------------------------------------------------------------------------------

   SUBROUTINE cheb_build_Fij(nnbox,nrho,npot,t,Fij)

      USE sefmm_quadrature, ONLY: get_roots_and_weights
      USE sefmm_cheb_utils, ONLY: cheb_get_nodes

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nnbox, nrho, npot
      REAL(REALK),   INTENT(IN)  :: t(:)
      REAL(REALK),   INTENT(OUT) :: Fij(:,:,:)

      REAL(REALK)   :: Vpts(npot)
      REAL(REALK)   :: rts(quad_pts), wts(quad_pts)
      REAL(REALK)   :: Tn(quad_pts,nrho)
      REAL(REALK)   :: T0, T1
      INTEGER(INTK) :: i,j, nt
      LOGICAL :: fail

      CALL CPU_TIME(T0)
      Fij(:,:,:) = zero

      ! Get potential grid points
      CALL cheb_get_V_nodes(nnbox,npot,Vpts)

      ! For numerical integration of F-vectors
      CALL get_roots_and_weights(quad_pts,quad_pts,-one,one,rts,wts)
      DO nt = 1, quad_pts
         CALL get_Tn(nrho,rts(nt),Tn(nt,:))
      END DO

      ! Build Fij for all t-points
      DO nt = 1, SIZE(t)
 !        IF ((nnbox > 0) .AND. (t(nt) > 100d0)) CYCLE  ! All v-points vanish
         IF (t(nt) > T_VALUE_THR) THEN
            CALL get_Fij_analytic(nrho,npot,t(nt),Vpts,Fij(:,:,nt))
         ELSE
            CALL get_Fij_numeric(nrho,npot,t(nt),Vpts,Tn,rts,wts,Fij(:,:,nt))
         END IF
      END DO

      CALL CPU_TIME(T1)
!      CALL TIMTXT('>>> TIME USED get F-vectors', T1-T0, LUPRI)

   END SUBROUTINE cheb_build_Fij

!------------------------------------------------------------------------------

   SUBROUTINE cheb_get_V_nodes(nnbox,npot,Vpts)

      USE sefmm_cheb_utils, ONLY: cheb_get_nodes

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nnbox, npot
      REAL(REALK),   INTENT(OUT) :: Vpts(:)

      INTEGER(INTK) :: tmp

      SELECT CASE (nnbox)
         CASE (0)
            CALL cheb_get_nodes(npot,Vpts)
         CASE (1)
            tmp = npot/3
            CALL cheb_get_nodes(tmp,Vpts(1:tmp))
            Vpts(tmp+1:2*tmp) = Vpts(1:tmp)
            Vpts(1:tmp) = Vpts(tmp+1:2*tmp) - two
            Vpts(2*tmp+1:3*tmp) = Vpts(tmp+1:2*tmp) + two
         CASE DEFAULT
            CALL fmm_quit('cheb_get_V_nodes')
      END SELECT

   END SUBROUTINE cheb_get_V_nodes

!------------------------------------------------------------------------------

   SUBROUTINE get_Fij_analytic(nrho,npot,t,Vpts,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nrho, npot
      REAL(REALK),   INTENT(IN)  :: t
      REAL(REALK),   INTENT(IN)  :: Vpts(npot)
      REAL(REALK),   INTENT(OUT) :: Fij(0:nrho-1,npot)

      REAL(REALK)   :: s
      INTEGER(INTK) :: k

      DO k = 1, npot       ! V-point
         s = Vpts(k)
         Fij(0,k) = get_F0(nrho,s,t)
         IF ( ABS(Fij(0,k)) < 1d-10 ) CYCLE  

         Fij(1,k) = get_F1(nrho,s,t)
         IF (nrho > 2) Fij(2,k) = get_F2(nrho,s,t)
         IF (nrho > 3) Fij(3,k) = get_F3(nrho,s,t)
         IF (nrho > 4) THEN
            CALL get_Fn(nrho,s,t,Fij(:,k))
         END IF

      END DO

   END SUBROUTINE get_Fij_analytic

!------------------------------------------------------------------------------

   SUBROUTINE get_Fij_numeric(nrho,npot,t,Vpts,Tn,rts,wts,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nrho, npot
      REAL(REALK),   INTENT(IN)  :: t
      REAL(REALK),   INTENT(IN)  :: Vpts(npot)
      REAL(REALK),   INTENT(IN)  :: Tn(:,:)
      REAL(REALK),   INTENT(IN)  :: rts(:), wts(:)
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      REAL(REALK)   :: s,fac
      INTEGER(INTK) :: n,k,nt

      Fij = 0d0
      DO k = 1, npot                ! V-point
         s = Vpts(k)

         DO n = 1, 1                ! Tn
            DO nt = 1, SIZE(rts)    ! GL-quadrature
               fac = t*(s-rts(nt))
               fac = -fac*fac
               Fij(n,k) = Fij(n,k) + wts(nt)*Tn(nt,n)*EXP(fac)
            END DO
         END DO

         IF ( ABS(Fij(1,k)) < 1d-10 ) CYCLE

         DO n = 2, nrho             ! Tn
            DO nt = 1, SIZE(rts)    ! GL-quadrature
               fac = t*(s-rts(nt))
               fac = -fac*fac
               Fij(n,k) = Fij(n,k) + wts(nt)*Tn(nt,n)*EXP(fac)
            END DO
         END DO

      END DO

   END SUBROUTINE get_Fij_numeric

!------------------------------------------------------------------------------

   SUBROUTINE get_Tn(npts,x,Tn)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: npts
      REAL(REALK),   INTENT(IN)  :: x
      REAL(REALK),   INTENT(OUT) :: Tn(0:npts-1)

      INTEGER(INTK) :: n

      Tn(0) = 1d0
      Tn(1) = x
      DO n = 1, npts-2
         Tn(n+1) = 2d0*x*Tn(n) - Tn(n-1)
      END DO

   END SUBROUTINE get_Tn

!------------------------------------------------------------------------------

   FUNCTION get_F0(npts,s,t)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: npts
      REAL(REALK),   INTENT(IN) :: s,t
      REAL(REALK) :: get_F0

      get_F0 = ROOTPI * (ERF(t*s +t) - ERF(t*s -t)) / (2d0*t)

   END FUNCTION get_F0

!------------------------------------------------------------------------------

   FUNCTION get_F1(npts,s,t)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: npts
      REAL(REALK),   INTENT(IN) :: s,t
      REAL(REALK) :: get_F1

      REAL(REALK) :: A,B1,B2, U,V

      A  = ROOTPI * (ERF(t*s +t) - ERF(t*s -t))
      B1 = EXP( -t*t*(s+1d0)*(s+1d0) )
      B2 = EXP( -t*t*(s-1d0)*(s-1d0) )

      U = (B1-B2) / (2d0*t*t)
      V = A * s/(2d0*t)

      get_F1 = U + V

   END FUNCTION get_F1

!------------------------------------------------------------------------------

   FUNCTION get_F2(npts,s,t)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: npts
      REAL(REALK),   INTENT(IN) :: s,t
      REAL(REALK) :: get_F2

      REAL(REALK) :: A,B1,B2, U,V

      A  = ROOTPI * (ERF(t*s +t) - ERF(t*s -t))
      B1 = EXP( -t*t*(s+1d0)*(s+1d0) )
      B2 = EXP( -t*t*(s-1d0)*(s-1d0) )

      U = ((s-1d0)*B1 - (s+1d0)*B2) / (t*t)
      V = A * (2d0*t*t*s*s + 1d0 -t*t) / (2d0*t**3)

      get_F2 = U + V

   END FUNCTION get_F2

!------------------------------------------------------------------------------

   FUNCTION get_F3(npts,s,t)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: npts
      REAL(REALK),   INTENT(IN) :: s,t
      REAL(REALK) :: get_F3

      REAL(REALK) :: A,B1,B2, U,V,U1,U2

      A  = ROOTPI * (ERF(t*s +t) - ERF(t*s -t))
      B1 = EXP( -t*t*(s+1d0)*(s+1d0) )
      B2 = EXP( -t*t*(s-1d0)*(s-1d0) )

      U1 = B1*( 4d0 -4d0*t*t*s +t*t +4d0*t*t*s*s)
      U2 = B2*(-4d0 -4d0*t*t*s -t*t -4d0*t*t*s*s)
      U  = (U1+U2) / (2d0*t**4)

      V = A * s*(4d0*t*t*s*s + 6d0 -3d0*t*t) / (2d0*t**3)

      get_F3 = U + V

   END FUNCTION get_F3

!------------------------------------------------------------------------------

   SUBROUTINE get_Fn(npts,s,t,Fn)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: npts
      REAL(REALK),   INTENT(IN)  :: s,t
      REAL(REALK),   INTENT(OUT) :: Fn(0:npts-1)

      REAL(REALK)   :: A,B,C,D,E,E1,E2
      INTEGER(INTK) :: n

      E1 = EXP( -t*t*(s+1d0)*(s+1d0) )
      E2 = EXP( -t*t*(s-1d0)*(s-1d0) )

      DO n = 2, npts-3

         A = (2d0*s*n -2d0*s)/(n-1d0)
         B = 2d0*( 1d0/(n-1d0) + (n*n-1d0)/(t*t*(n-1d0)) )
         C = 2d0*(-s*n-s)/(n-1d0)
         D = (n+1d0)/(n-1d0)

         E = (2d0/(t*t*(n-1d0))) * (E1*(-1d0)**n +E2)
     
!print '(I4,5F14.6)', n+2,A,B,C,D,E
         Fn(n+2) = A*Fn(n+1) +B*Fn(n) +C*Fn(n-1) +D*Fn(n-2) +E

      END DO
!do n = 0, npts-1
!print *, n, Fn(n)
!end do

   END SUBROUTINE get_Fn

!------------------------------------------------------------------------------

END MODULE sefmm_cheb_Fij
